In [2]:
%pylab inline
import random
The method of metadynamics seeks to explore a free energy landscape and cross high energy barriers using artificial bias potentials. This allows a system to pass through a highly unfavorable state without needing to quench the system. This is done by adding Gaussian potentials in deep energy wells, filling them, and allowing the system to "spill" into another metastable state.
$$ V_{bias} = G(s)e^{-\frac{(s-s_0)^2}{2\sigma ^2}} $$
In [3]:
figure(figsize=[6,2])
height = 5
center = 0
width = 2
x = linspace(-10,10,1000)
gaussian = height*exp(-(x-center)**2/(2*width)**2)
plot(x,gaussian)
title("Gaussian")
Out[3]:
As the system moves forward in time, a Gaussian is added to the free energy surface $F(s)$ at the previous location.
$$ F'(r) = F(r) + \sum G(s)e^{-\frac{(s-s_0)^2}{2\sigma ^2}} $$The plot below shows the effect of this bias potential on the surface. As more gaussians are added to the system, wells are filled and the MD simulation is able to sample the whole system.
In [3]:
figure(figsize=[10,5])
x = linspace(-1,1,1000)
height = .2
center = 0
width = .02
gaussian = height*exp(-(x-center)**2/(2*width)**2)
y = (3*x**10)+(2*x**9)-(x**8)+(3*x**7)-(2*x**6)-(5*x**5)-(6*x**4)+(x**3)+(7*x**2)-(x)
y2 = (3*x**10)+(2*x**9)-(x**8)+(3*x**7)-(2*x**6)-(5*x**5)-(6*x**4)+(x**3)+(7*x**2)-(x) + gaussian
fig = plt.figure(1)
plot(x,y, 'r')
plot(x,y2,'b')
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')
Out[3]:
This cell demostrates the behavior of an MD simulation that uses Metadynamics to sample all of surface. Notice there are a few key parameters that determine the computational expense and resolution of the results:
Depending on how rough the free energy landscape is, these parameters will need to be adjusted to smooth the surface and sample key metastable states.
In the example below, the red line describes the "real" free energy landscape (same arbitrary function for this example). The blue line is the sum of the bias potentials added to the free energy landscape. All the energy wells are sampled in this system and "filled" with artificial potentials.
In [4]:
figure(figsize=[10,6])
def y(x):
return (3*x**10)+(2*x**9)-(x**8)+(3*x**7)-(2*x**6)-(5*x**5)-(6*x**4)+(x**3)+(7*x**2)-(x)
def add_gaussian(height, width, func):
gaussian = height*exp(-(x-center)**2/(2*width)**2)
return func(x) + gaussian
x = linspace(-1,1,1000)
surface = y(x)
y1 = surface
plot(x,surface, 'r')
height = .06
center = 0
width = .01
Gaussian = 0
for i in range(0,30000):
center = random.uniform((-1*.2)+center,(1*.2)+center)
gaussian = height*exp(-(x-center)**2/(2*width)**2)
y2 = y1 + gaussian
if max(y2) > 3:
pass
else:
y1 += gaussian
Gaussian += gaussian
plot(x,y1, 'b')
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')
Out[4]:
Take the total artificial bias potential (sum of the Gaussians) and plot it's negative
$$ F(r) \approx - \sum G(s)e^{-\frac{(s-s_0)^2}{2\sigma ^2}} $$Depending on how robust and resolved the metadynamics simulation is, the free energy landscape is approximately equal to this bias potential.
In [22]:
figure(figsize=[12,8])
plot(x,-Gaussian+3.1)
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')
figure()
figure(figsize=[12,8])
plot(x,-Gaussian+3.1)
plot(x,y(x), 'r')
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')
Out[22]:
The importance of these type of free energy landscapes is highly dependent on the choice of collective variable's (a.k.a. reaction coordinates). There are three major criteria for CV's:
The collective variables are the key to determining the mechanism for transitions between metastable states. Metadynamics takes a "course-grained" approach to sampling a system and finding these metastable states, focusing only on "interesting" variables that describe the dynamics of the system. These variables need to be chosen correctly to find results with meaning.